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Abstract 

We describe a Lohner-type algorithm for the computation of rigorous 
upper bounds for reachable set for control systems, solutions of ordinary 
differential inclusions and perturbations of ODEs. 

1 Introduction 

Our goal is to present a Lohner-type algorithm for an rigorous integration of 
perturbations of ODEs, which can be seen also as an algorithm for an integration 
of control systems or ordinary differential inclusions. This paper depends heavily 
on [Zl], as the proposed algorithm is a modification running on top of the C°- 
Lohner algorithm for ODEs described (after |Lo[ ILolj ) there. 
We study the following nonautonomous ODE 

x'(t) = f(x(t),y(t)), x(0)=x o (1) 

where x G R™, / : R™ x R m ->■ R n is C 1 and y : R D D -> W n . Assume that 
we have some knowledge about y(t), for example \y(t)\ < e for < t < T. We 
would like to find an rigorous enclosure for x(t). 

The problem of this type arises, for example, in the context of the control 
theory (see [GlIKSllSz] ) and in the rigorous integration of dissipative PDEs (see 
|ZM[ IZ21 IZ4] for more details). In this last setting x represents the dominating 
modes and y is a tail of the Fourier expansion, so that ([T]) is complemented 
by the equation for y of the form y'(t) = g(x(t), y(t)) for which we are able to 
produce some a priori bounds. The proposed algorithm works, as we were able 
using it prove the existence of multiple periodic orbits for Kuramoto-Sivashinsky 
PDE [Z21IZ4] . 

The proposed algorithm can also be used to find rigorous bounds for solutions 
of differential inclusions 

x' e h{x) + e(t), (2) 

where h : R™ -> R" is a C^vector field and e(t) C R". We can cast (0) in the 
form ([T]) by setting f(x, y) — h(x) + y and requiring that y(t) € e(t) for all t. 
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Non-autonomous ODEs represent another important class of applications. 
While one can easily modify the Lohncr algorithm to handle a non-autonomous 
ODE directly, it makes sense to apply the proposed Lohner-type algorithm for 
perturbed ODEs for (fTJ) , because only in this way we can estimate rigorously the 
Poincare map on a section ot(x) = (defined in terms of x only) for any initial 
conditions (a:, to)- This kind of algorithm shall allow to attack the question of 
symbolic dynamics for non- autonomous ODEs (see |CZ] ) and ODEs with small 
delays (see |WjZ| ). 

Another new element in this paper, besides the proposed algorithm, is a 
new inequality concerning bounds for perturbations of ODEs. It is contained in 
Theorem Inland links together the component- wise estimates based on one-sided 
Lipschitz conditions (see [W]) and the logarithmic norms (see [Dl |L]V 

The content of the present paper can be described as follows: in Section[2]we 
define a notion of weak solution of (fl]) and state some facts from the theory of 
Lebesgue integration. In Section [3] we recall the notion of the logarithmic norm 
and state its basic properties. In Sections 0] and [5] we derive basic estimates 
for comparison of perturbed and unperturbed ODEs. In Section [5] we give a 
description of one step of the proposed Lohner-type algorithm. In Section [7] we 
describe how to estimate the trajectory of |T]) between time steps which allows 
to compute the Poincare map. In the following section we discuss some tests. 

The algorithm presented in this paper was implemented as a part of CAPD 
library (see [CAPD] ). This library contains many tools for rigorous computa- 
tions and computer assisted proofs in the contexts of dynamical systems. All 
the tests in Section [8] was performed using CAPD library. 

1.1 Basic notation 

We will use the same conventions as in |Zlj . In the sequel, by arabic letters 
we denote single valued objects like vectors, real numbers, matrices. Quite 
often in this paper we will use square brackets, for example [r], to denote sets. 
Usually this will be some set constructed in the algorithm. Sets will also be 
denoted by single letters, for example S, when it is clear from the context that 
it represents a set. In situations when we want to stress (for example in the 
detailed description of algorithm) that we have a set in a formula involving both 
single-valued objects and sets we will rather use the square bracket, hence we 
prefer to write [S] instead of S to represent a set. From this point of view [S] 
and S are different symbols in the alphabet used to name variables and formally 
speaking there is no relation between the set represented by [S] and the object 
represented by S. Quite often in the description of the algorithm we will have 
a situation that both variables [S] and S are used simultaneously, then usually 
S € [S], but this is always stated explicitly. 

For a set [S] by [S]i we denote the interval hull of [S], i.e. the smallest 
product of intervals containing [S]. The symbol hull(xi, . . . , Xk) will denote the 
interval hull of intervals x\, . . . , Xk- For any interval set [S] = [S]i by m([5 1 ]) we 
will denote a center point of [S]i. For any interval [a, b] we define a diameter 
by diam([a, b]) — b — a. For an interval vector or an interval matrix [S] = [S]i 
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by diam ([S]) we will denote the maximum of diameters of its components. For 
an interval [£ _ ,2H-] we set right([x~ , x + ]) = x + and left([x~, x + }) = x~. 

If f(x%, X2, ■ ■ ■ , Xk) is a function and let X\, X2, ■ ■ . , Xk be some sets, then 

by 

f(X%, . . .,X k ) = {f(zx, ...,Zk)\ where z l £ X l for % = 1, . . ., j} 

For a set X c R d by int X we denote an interior of X . For R™ we will denote 
the norm of x by ||x|| and if the formula for the norm is not specified in some 
context, then it means that it is ok to use any norm there. Let xq £ K s , then 
B(x Ql r)={zeR s I Ha* -z\\ <r}. 

For v,w eM. n and A, B G R nxn (n = 1, . . . , 00) we say that 

v < w iff Vi Vi < Wi, 

A < B iff Vij AijKBij. 

1.2 Warning. 

At the first encounter with the question of an rigorous integration of ([T]) one 
may hope that the direct application of any algorithm for rigorous integration 
of ODEs should be enough for @. To this end consider a differential inclusion 

x'ef(x) + [e], [e] = IIg =1 [-e i ,e i ]. (3) 

and a related ODE 

x' = f(x)+e, ee[e]. (4) 

One may naively hope that, for example, the Lohner algorithm applied to 
l[4"j) with [e] as an interval parameter in the definition of a constant term in f{x) 
will give an enclosure not only for ([4]), but also for ([3]). For this to be true we 
need the following 

Conjecture 1 Assume x{t) satisfies for t £ [0, T] . 

Then for any t S [0, T] there exists e £ [e] such that x e (t) = x(t) and 
x e (0) = x(0), where x e is a solution of 0). 

The above conjecture is false as shown by the following example [Sej . 
Consider a differential inclusion given by 

x' £ y+[-e,e], (5) 
y' € -x+[-e,e\. 

For fixed 8 £ [— e, e] 2 we have the following system of ODEs 

x' = y + 81, (6) 
y' = -x + S 2 , 

all solutions with an initial condition in a compact set have a uniform bound 
independent of 6 for t > 0, which is given by the energy integral for (|6|) 

(x-S^ + iy + S^ 2 . (7) 
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This is not the case for the solutions of ([5]) as it is clearly seen for e(t) given as 
a resonant forcing 



x 



-x + esini. 



(8) 



2 Control Systems, the notion of the solution 

In this section we define a notion of (weak) solution of ([T]). 

We use some standard notions from the measure theory, see for example [Ruj 
for precise definitions. The integral will always mean the Lebesgue integral and 
the measure of the set is always Lebesgue measure. 



2.1 Some facts from the theory of Lebesgue integral 

We will denote by m(E) the Lebesgue measure of E. 

Let D be a measurable subset of R fe . By L 1 (D) we will denote a set of 
measurable functions / : D — > R such that J D \f\dm < oo. If / : D — > R" is 
measurable, then we say that / e L 1 (D) if function ||/|| £ L 1 (D). 

Definition 1 Let D C R be an interval. Function f : D — > R fc is absolutely 
continuous , i/ for every e > there exists S > 0, such that for any family of 
disjoint intervals (ai,/3i), . . . , (aAr,/3jv) such that 

N 
i=l 

i/ie following inequality is satisfied 

N 

EW)-/(«o)<e 

i=l 

The following statement follows directly from results about the differentia- 
bility of measures and functions of bounded variation (see |Ru[ Chapter 8]). 

Theorem 2 Let D = [a, b], x ; D -> R n . 

There exists g : D — * R" a measurable function such that equation 

x{t) - x(a) = I g{s)ds (9) 

/io/rfs /or aZ/ i £ [a, 6] i/f x is absolutely continuous. In this situation x'(t) exists 
almost everywhere in [a, b] and x'(t) = g(t). 
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Definition 2 Assume x £ M. k . We say that a sequence {E{\ of measurable 
subsets of M. k converges well to the point x, if there exists a > such that: 
every set Ei is contained in B(x,Vi), such that 

m{Ei) > am{B[x,ri)), lim r.; = (10) 

i — >oo 

In the sequel we will need the following theorem 

Theorem 3 JRul Thm. 8.8] Assume that f £ L 1 ^) and defi ne a Lebesgue 
set Lf of the function f as the set of all points xq £ M. k for which 

lim — / |/(a;) - f(x )\dx = (11) 

for every sequence {Ei} converging well to the point Xq. 
Then set Lf contains almost all points ofR k . 

The above theorem immediately implies the following lemma. 

Lemma 4 Let f : [a, b] — > M fc be a measurable function. Then for almost all 
points x £ [a, b) holds 

Km \ r +h \\f( S )-f(x)\\ds = (12) 

2.2 Weak solutions of ODEs 

Control System is given by equation 

x'{t)=f(x(t) >y (t)) x(t ) = x (13) 

where x £ E n , / : R" X W n -> M" is C 1 and j/:M3J}^ W n is a measurable 
function from a given class U . 

Because the right hand side of (|13p can be non-continuous we need to define 
what we mean by solution of (fT3f . 

Definition 3 Let D C K be an interval (a connected subset ofM.) containing 
to- 

An absolutely continuous function x : D — ► R n is a weak solution of \13\) if 
for all t £ D holds 

x(t)=x + [ f(x{s),y(s))ds. (14) 

J t 

We say that a continuous function x : D — * K.™ is a (classical) solution of 
H13\) if x'(t) exists for all t £ miD, x(to) — xq and 

x'(t)=f(x(t),y(t)), Vt£intD. (15) 
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From Theorem [5] it follows that a: is a weak solution of (fT5)) iff 

x'(t) — f(x(t),y(t)), allmost everywhere in D (16) 

and the function t i— > f(x(t), y(t)) is in L}(D). Hence the weak solution in the 
sense of Def . [3] is a solution of (|13p in the sense of Caratheodory [W] . 

In the remainder of this paper we will always consider the function / on the 
right hand side of (fT5| to be of class C k (for k > 1) and y to be bounded on 
compact intervals and measurable. In such situation the integral equation (| 14[) 
has a unique solution defined for some h > on [to, t + h] . The proof of this 
fact is a straightforward application of the Banach contraction principle [W] . 



3 Basic facts on logarithmic norms 

Let || • || denote a vector norm on R™ as well as its subordinate matrix (operator) 
norm on R™ xn . The classical definition of the logarithmic norm of matrix A, 

\\I + hA\\-l 

= A — h — (17) 

was introduced in 1958 independently by Dahlquist [D] and Lozinskii [L]. 

In this section we will briefly recall some basic facts, with proofs, about the 
logarithmic norms. For survey regarding the modern developments stemming 
from this notion the reader is referred to [So] and the literature given there. 
Our presentation is based on |DV| Ch. 1.5 ], which was based on [D] . 

Lemma 5 For any matrix A S R™ x ". The limit in exists and 

+ < + for0< hl <k 2 (18) 

hi /12 

-Ml < v(A)<\\A\\. (19) 
Proof: Let us fix h > and let < 9 < 1, then 

\\I + 6hA\\ = \\9(I + hA) + (1 - 9)I\\ < 9\\I + hA\\ + (1 - 6»)||/||. 
From this immediately obtain 

\\I + 9hA\\-l < \\I + hA\\-l (2Q) 



9h h 

which proves l|18p. 

From the triangle inequality one gets 

-h\\A\\ < \\I + hA\\ - \\I\\ < h\\A\\, (21) 

therefore 

„ „ ||7 + hA\\ - 1 , s 

-Il4l<- <\\M- (22) 

The monotonicity (|18[) and the existence of the lower bound imply the existence 
of fi(A). I 
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Theorem 6 The function \i : R" xn — > R, which assigns to A its logarithmic 
norm is continuous and convex. Moreover, functions fJ,(h, A) = H^+^H-i CQn _ 
verge locally uniformly and monotonically to (J,(A) for h — > + . 

To be more precise, for any compact set K C R" xn and any e > there 
exists ho > 0, such that for all < h < ho and any A £ K holds 

e > fj,(h, A) - fi(A) > 0. (23) 

Proof: Let h > 0. An easy computation show that, for any < A < 1 and 
Ai,A 2 e R nx " holds 

ft(h, XAi + (1 - X)A 2 ) < Xfi(h, Ai) + (1 - X)n(h, A 2 ). 

Therefore, for any h > function fj,(h, •) : R" xn — > R is convex. 

By taking the limit h — > + from Lemma [5] it follows that is a con- 

vex function. Observe that on any bounded set U C R™ x ™ fi(A) is bounded 
by sup^gfj 1 1 1 < +oo, therefore from the theory of convex functions (see for 
example [Lai Chap. 6]) it follows that [i is continuous. The uniform conver- 
gence of fi(h, •) to /i on compact sets follows from Dini's Theorem on monotone 
sequences of pointwise converging continuous functions to continuous limit and 
Lemma [5] I 

The following lemma follows directly from the convexity of l-i(A) 
Lemma 7 Let A : [0, 1] — > R nx ™ be a bounded measurable function. Then 

A{s)ds\ < [ n(A{s))ds < sup (i,(A(s)). (24) 

\Jo ) JO s6[0,l] 



4 Bounds for perturbations of ODEs 

In this section we state the basic theorem comparing a solution of an ODE 
and an approximate solution. Our approach unifies the approach based on 
logarithmic norms and one-sided Lipschitz condition leading to component-wise 
bounds from [Wl Ch. 11.13]. 

4.1 Estimates for non-autonomous linear equations 

Consider a linear equation 

x'{t) = A(t) ■ x{t) + b{t), (25) 

where x(t) € R fc , A(t) € R fcxfc , b(t) g R fe , A and b are bounded and measurable. 

We would like give some bounds on solutions of (|25[) . We assume that our 
phase space R fc is decomposed as follows R fc = ©™ =1 R fci . Therefore, we have a 
decomposition of z € R fc into (z\, . . . , z n ) such that £ R fci . In this section we 
will carefully distinguish between the symbol |j • || and | • |. The symbol || • || will 
always denote a norm, but the symbol \z\ for z £ R fe will usually denote a vector 
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of norms of Z{, but this will be always clearly indicated in the text. Observe 
that, when we have such decomposition, then equation (|25|) can be written as 
follows 

z' i (t) = ^A ij (t)z j (t) + b i (t), i = l,...,n (26) 

3 

where Zuh <E R fe * and A i3 (t) e L(R fe *,R fe i) is a linear map (a matrix). In this 
way matrix A is decomposed into blocks Ay. For each block we will assign 
number Jy and collect them in matrix J. Roughly speaking Jy will estimate 
an influence of Zj on z[. 

The fundamental lemma in this section is: 

Lemma 8 Assume that z : [0,T] — > R fe = ©™ =1 R' Ci is an absolutely continuous 
map, which is a weak solution of the equation 

z'(t) = A(t)-z{t) + 5(t), (27) 

where S : [0, T] -> R k and A : [0, T] -> R fexfe are bounded and measurable. 

Assume that measurable matrix function J : [0,T] — > R™ xn satisfies the 
following inequalities for all t € [0, T] 

, . f II ^4i, (i) II for i ^ j, 
Jijit) > { (28) 
yn{A u {t)) fori = j. 

Let d(t) = \\5i(t)\\ and \z\(t) = {\\zi(t)\\, ||* 2 (*)||, • ■ ■ , IM*)||)- 
Then 

\z\(t)<y(t) (29) 
where y : [0, T] — > R" is a weafc solution of the problem 

y'(t) = J{t)y(t) + C(t), y(0) = |*|(0). (30) 

Proof: Observe that for all i the function t i— > ||zj(t)|| is absolutely continu- 
ous. Therefore from Theorem [5] it follows that for almost every t G [0, T] the 
derivative of \\zi\\ exists. We will estimate this derivative for such t. 
We have 

/t+h pt+h 
A(s)z(s)ds + / 6(s)ds = 

ft+h 

z{t) + h (A(t)z(t)) + hS(t)) + J (A(s)z(s) - A{t)z{t)) + (6(s) - S(t))ds 

Let us fix i and t £ [0,T). We consider the projection onto i-th subspace. 
We have 

\\zi(t + h)\\ < \\I + hA u (t)\\ ■ \\z t \\(t) + hJ2\\ A ij(t)\\ ■ \\xj(t)\\ + h\\8i(t)\\ + 

rt+h pt+h 

\\A(s)z{s)-A(t)z(t)\\ds+ / \\5(s)-S(t)\\ds 
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and then we obtain for h > 



||*i(* + /Q||-||*i(t)|| < \\I + hA u (t)\\~l + 
h ~ h 1 

J2 \\Mt)\\ ■ IkiWII + Ci + ~ P( s )z(,s) - A(t)z(t)\\ ds + 



t+h 

||<S(s)-*(t)||da 



Observe that from Lemma 0] it follows that the last two terms in the above 
inequality tend to as h — > for almost all points in [0, T). From now on we 
assume that t is such point. 

By passing to the limit with h — > + we obtain for almost all points in 

te[0,T] 

%%) < MM*))IWI(*) + E W A ^)\\ ■ \\zj(t)\\ + ci(t) < 



dt 



E^'WPillW + ^W ( 31 ) 



Let us define 

= (a:i(t),a:2(t) s ... ,«»(*)) = (Pi(i)ll, ll^(*)||, . • • , ||^(*)||), 

Inequality (|3 1 [) can be rewritten in vector form as follows 

x'(t) < J(t) ■ x(t) + C(t), for almost all t G [0, T]. (32) 

Let y : [0, T] — > R™ be a weak solution of 

y'(t) = J(t)-y(t)+C{t), (33) 

such that ?/(0) > |z|(0) = x(0). 
We want to show that 

x(t)<y(t), te[0,T]. (34) 

Let us take diagonal matrix A G R" xn , such that A„ + Ju(t) > for all i = 
l,...,n and t G [0, T]. Let us define matrix-valued function B : [0, T] -> R" x " 

by 

B{t)=k + J(t). (35) 

Obviously > for all f G [0, T]. 

For any i = 1, . . . , n from (f3"2"| we obtain for almost all t £ [0, T] 

^(t)+A ii a; i (t) <J2 B ij( t ) x j( t )+ C i( t )' ( 36 ) 
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hence 

| (e A «^(i)) < e A ^ fa Bax^t) + d(t) 

The last inequality has the following vector form 

i {e M x(t)) < e At B(t)x(t) + e At C(t). (37) 

From the above inequality and from Theorem [2] it follows that 

/"* d 

e At x(t) = e A -°x(0) + I j (e At x(t)) (s)ds < x(0) + 
ct 

e As B(s)x{s) + e As C(s)ds. 

Jo 

Hence we obtain 

x(t) < e- At x(0) + f e - Ht - s) (B(s)x(s) + C(s)) ds for t € [0,T] (38) 



o 



An analogous computation applied to p3|) shows that y satisfies the following 
integral equation 

y(t) = e- At y(0) + f e~ A ^ (B(s)y(s) + C(s)) ds. (39) 



o 



Now we are ready to prove (|34|) . Let 

to = sup{t e [0, T] I y(s) > x(s), s € [0, «)}. (40) 

Obviously from the continuity of y(t) — x(t) it follows that to > 0. From (|39p 
and (|3"81 we obtain 



(to) - x(t ) > e~ At °(y(0) - x(0)) + / e- A ^-^B(s)(y(s) - x(s)) ds > 0. 



y( 

By the continuity inequality y(t) > x(t) will hold for t € [to, to + e) for some 
e > 0. Therefore t = T. 

Hence condition (f34|) holds. By passing to the limit y(0) — > x(0) we obtain 
our assertion. I 

Theorem 9 Let ft, > 0. Assume that f : R n x M m -> R" 6e C 1 and j/ : 

[to, to + ft] — » K m *s bounded and measurable. 

Let [W y ] C W n be convex and such that, y([to,to + ft]) C [W y ]. 
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Let y c € [Wy\. Assume that x\,x 2 ■ [to, to + h] - * ^ n , both absolutely 
continuous, are weak solutions of the following problems, respectively 

x 'i = f(xi,y c ), Xi(t )=x , (41) 
x' 2 = f(x2,y(t)), x 2 (t )=x . (42) 

Let \W\] C [W 2 ] C M. n be convex and compact and such that 

Xi(t)e[Wx), x 2 (t)e[W 2 ], fort € [to, to + ft]- 

Then the following inequality holds for t G [to, to + h] and i = 1, . . . , n 

\x 1 , i (t)-x 2 , i (t)\ < (^.(.o-xo^+j^^Cdsj , 



(43) 



where 



[S] = {f(x ) y c )-f(x,y)\xe[W 1 ],ye[W y ]}, 
Ci > sup \[6i]\ , i = l,...,n 



Jij > 



supt,m([W 2 },[W y ])) if ■i = j 



dxj 

ifi + 3- 



sup 

Proof: Let z(t) = xi(t) — x 2 {t). We have for t G [to, to + h] 



z(t) = ( xi(*o) + J f{xi(s),y c )ds J - yx 2 (t ) + J f(x 2 (s),y(s))ds 



z ( t o) + / (f{xi(s),y c )-f(x 2 {s),y(s)))ds. 



to 



Now observe that 



f{x x {t),y c ) - f(x 2 (t),y(t)) = f( Xl (t),y c ) - f( Xl (t), y(t)) + 
/(xi (<),»(*)) - f(x 2 (t),y(t)) = + A(t) ■ (a:i(*) - x 2 (t)), 

where S(t) € [5] is bounded and measurable and 

Aij{t) = jf ^-(ar 2 (t)+a(a;i(t)-ar2(t)),y(t))da 

is bounded and measurable matrix. 
We obtain ^ 

z(t)=z(t )+ [ (A(s)z(s)+S(s))ds (44) 

To apply Lemma [8] to the function z = xi — x 2 to obtain (|51[) we need to 
show that 

j., > | SU P*G[to,to+ft] II^WII fOT l + 3, (45) 
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For the off-diagonal terms we have 



IIAiWII < 



1 df 

(x 2 (t) + s( Xl (t) - x 2 (t)),y(t)) 



sup 

xE[W 2 ],ye[W v ] 



dfi 
dxi 



ds < 



< Ji, 



For the diagonal case we use Lemma [7J 

The result now follows from Lemma [5] I 

It is possible to organize the error estimates slightly differently, namely esti- 
mate [S] on [W2] x [W v ] instead of on [Wi] x [W y ], which will produce larger [5], 
but in the same time estimate J on [W2] x {y c } instead of [W2] x [W y ], which 
should result in better J, to obtain the following variant of the above theorem. 

Theorem 10 The same assumptions and notations as in Theorem^ 
Then the following inequality holds for t E [to, to + h] and i = 1, . . . , n 

\xxM -x 2 ,i{t)\ < (e^-^ ■ Ozo-zo)). + e J{t ~ s) Cds^j , 



(46) 



where 



[5] = {f(x,y B ) - f(x,y) \ x G [W 2 ],y G [W v ]}, 
Ci > sup I [<5i] I , i = l,...,n 



J j 



> 



Bupn(§£;(\W 2 ],y e )) ifi = j, 



sup 



dx-j 



ifi + J- 



Proof: We proceed as in the proof of Theorem [9] But the difference between 
f(xi(t),y c ) and f(x2(t),y(t)) is computed differently. Namely, 

f{x x {t),y c ) - f{x 2 (t),y{t)) = f( Xl {t),y c ) - f(x 2 (t),y c ) + 
f(x 2 (t),y c ) - f(x 2 (t),y(t)) = A(t) ■ {x x {t) - x 2 {t)) + S(t), 



where S(t) G [5] and 



1 df, 



A y( < )= / ir ± (Mt) + s(x 1 {t)-x2{t)),y c )ds. 



1 dxj 

We continue as in the proof of Theorem [9] 



5 Formulas for various cases 

In this section we rewrite Theorems [9] and [10] in the form, which will be later 
used in our algorithm for the integration of differential inclusions. 
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5.1 The estimation of perturbations of ODEs based on 
logarithmic norms 

From Theorem 1101 using the trivial decomposition consisting of the whole space 
we obtain the following lemma. 

Lemma 11 Let h > 0. Assume that / : R™ x R m — > R" be C 1 and y : 
[to, to + h] — > R m be bounded and measurable. 

Let [W y ] C R m be convex and such that, y([io,*o + h]) C [W y ]. 

Let y c £ [W y ] . Assume that Xx,x 2 ■ [to, to + h] — > 1" both absolutely contin- 
uous, are weak solutions of the following problems, respectively 

x'i = f(xi,y c ), Xx(to)=x , (47) 
x' 2 = f{x 2 ,y(t)), x 2 (t )=x . (48) 

Let [Wi] C [W 2 ] C 1" be convex and compact and such that 

x 1 (t)€[W 1 ], x 2 (t)e[W 2 ], fors£[t ,to + h]. 

Then for any t € [0, h] holds 

\\x 2 (t + t) - Xl (t + t)\\ < 

rto+t 

exp(lt)\\xi(t ) - x 2 (t )\\ + exp(lt) / exp(— Zs) ||[<5]||ds = 

J to 

exp(It)||asi(to) - aj 2 (*o)|| + M(exp(Zt) - 1) 

where I = sup ^(§£([^2], 2/c))J > M * s ^ e logarithmic norm of the matrix 
(see [HNWl for the definition) and 

[S] = {f(x,y c ) - f(x,y) \ x & [W 2 ],y e [W y ]}. 

5.2 A component- wise estimate 

From Theorem [§] using the trivial decomposition R m = we obtain the 

following lemma. 

Lemma 12 Let h > 0. Assume that / : R™ x R m — > R" be C 1 and y : 
[to, to + h] — > R m is bounded and measurable. 

Let [W y ] C R m &e convex and such that, y([to,to + h]) C [Wj,]. 

Lei y c € [Wy]- Assume that X\,x 2 : [to, to + h] — > R™, &oi/i absolutely 
continuous, are weak solutions of the following problems, respectively 

x'i = f{xi,y c ), xi(t )=x , (49) 
z 2 = f(x2,y(t)), x 2 (t ) = x . (50) 

Let [Wi] C [W2] C R" be convex and compact and such that 
ii(t)e[Wi], i2(t)e[M^], /oraG [io.to + Zi]. 
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wher 



Then the following inequality holds for t £ [to, to + h] and i = 1, . . . , n 

\x hi (t) - sc 2 ,*(*)] < (e Jt ■ (x - x )). + Qf e J C'-">C ds) , (51) 

[<5] = {/(^yc)-/(^,y)ke [W 1 ],y£ [W y ]}, 

Ci > sup | [<5i] | , i = l,...,n 

sup^a^a],^]) ifi = j, 



Jij > 



sup 



6 The Lohner-type algorithm for perturbations 
of ODEs 

For a given measurable and bounded on compact intervals function y : [0, oo) — > 
R m let ip(t, Xo, y) denotes a weak solution of equation (fl]) with initial condition 
x(Q) = x . For a given y £ R m let Tp(t,x ,y ) be a solution of the following 
Cauchy problem 

x' = f(x,y ), x(0)=x o (52) 

with the same initial condition a;(0) = xq. Observe that system (|52[) is a par- 
ticular case of dJ) with y(t) = y - 

Let U be a some family of functions y : [0, oo) — » M™ 1 which are measurable 
and are uniformly bounded on any compact interval, i.e. for any T > there 
exists M (T), such that for every y £ U and every t £ [0,T] holds ||y(i)|| < M(T). 

We are interested in finding rigorous bounds for cj)(t, [xq], [yo]), where [xq] C 
R™ and [yo] C [/. The set [yo] might be defined as some dynamical process, in 
this case we may need to compute something for each time step, or it can be 
just given by the specifying the bounds, for example y £ [yo] iff y(i) £ [— e, e] m 
and y is measurable. 

Below we propose a modification of the original Lohner algorithm |Lo| ILolj 
to treat problem (flj . Our presentation follows the description of the C°-Lohner 
algorithm presented in |Z1| . 



6.1 One step of the algorithm 

In the description below the objects with an index k refer to the current values 
and those with an index k + 1 are the values after the next time step. 
We define 

[Vk] = {y £ U | y(t) = z(t k + 1) for some z £ [y Q ]}. 
For given [y] C U we will also use the following notation 

[y]{[t 1 ,t 2 ]) = {z{t)]z£[y],t£[t 1 ,t 2 ]}. 
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One step of the Lohner algorithm is a shift along the trajectory of system 
{T]) with following input and output data: 
Input data: 

• tk is a current time 

• hk is a time step 

• [xk] C M™, such that ip(t k , [x ], [yo]) C [a;*] 

• eventually some bounds for [y^] 
Output data: 

• tk+i = tk + hk is a new current time 

• [xfc+i] C M", such that v(f fc+1 , [x ], [y ]) C [x fc+1 ] 

• eventually some bounds for [yo] [0, tk+i)- 

We do not specify here a form (a representation) of sets [x&]. They can be 
interval sets, balls, doubletons etc. (see |MZ| IZlj ), This issue is very important 
in handling of the wrapping effect and is discussed in detail in |Lo|, ILolj (see 
also Section 3 in |Zlj). 

One step of the algorithm consists from the following parts: 

1. Generation of a priori bounds for ip and [j/o] ([^fej tfe+l])- 

We find a convex and compact set [W2] C M. n and a convex set [W y ] C M m , 



2. We fix y c g [W y ], 

3. Computation of an unperturbed x-projection. We apply one step 

of the C°-Lohner algorithm to ([52| with a time step hk and an initial 
condition given by [xk] and yo — y c - As a result we obtain [xk+i] C R n 
and a convex and compact set [WjJ C K n , such that 



4. Computation of the influence of the perturbation. Using formulas 
from Lemmas IT21 or [TT| we find a set [A] C K n , such that 



such that 



<p([o,h k ]Mhk])c\w 2 ] 

[y k }([0,h k })G[W y } 



(53) 
(54) 



<p{hk,[xk],y c ) c [x k +i] 
V([0,h k ],[xk],Vc) c [W^i] 



<p(tk+i, [xo], [yo]) C <p(/i fe , [xk],y c ) + [A]. 



(55) 



Hence 



y(ifc+i, [x ], [yo]) C [xai+i] = [xk+i] + [A] 



(56) 



5. Eventually we do some computation to obtain [y/c+i] 
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6.2 Part 1 - comments 



In the context of an nonautonomous ODE with small and uniformly bounded 
[S] we can set [W y ] = M. To obtain [W 2 ] any rough enclosure procedure devised 
for ODEs should work. In the context of a dissipative PDE the whole story is 
more complicated and we refer the interested reader to |Z2j . 

6.3 Part 4 - details 

In Lemmas [TT1 and [T2l we have presented two ways to compute [A] = [A](/i) for 
< h < h k . 

An approach based on component-wise estimates 



1. We set 



[8] = [{f(x,y c )-f(x,y)\x&[W 1 ],ye[W y ]}] I 
d = right(|[£i]|), i = l,...,n 

gU(p-{[W 2 i\W v })) Hi=j, 

right (j§§-([W 2 ],[W y ])\). ifi^j. 



2. D = ^e Ji - h -^Cds 

3. [A,;] = [-Di, Di], for i = l,...,n 

It remains to explain how we compute f Q e A ^ s ^C ds. First observe that 



(At) 

[n + ' 

We fix any norm || • ||, such that for any matrix A = (a.y) we have |<Zjj| < \\A\\. 
It is not true for general norm, for example if we take vector norm on R 2 
defined by |j (x\ , x 2 ) | = max{ j^q^i , x 2 } then associated matrix norm of a matrix 

is equal to 1. We take for example L°°-norm, i.e. ||a;||oo = max^ | ar^ | 





(we should rather chose a norm for which \\At\\ is the smallest one). Let us set 

A m 

A — At^ Ay 



(m+1)! 
In this notation 



~ {A tr ^ , 

^ (m + l)! ^ " 

m=0 v ; m=0 



A Q = Id, A rn+ i = Am ■ — 

m + 2 



1G 



For the remainder term we will use the following estimate 

k 

I 

\A N+k \\ < \\A N \ 



N + 2 



Hence if 



N+2 



< 1, then 

E A 



m>N 



< 



\A 



N 



A 



= \\A 



N 



N + 2 



1 - 



.4 



N + 2 



And finally, 



N + 2- \\A\\ 



N 



Y^A m + [-r,rr (58) 

rn— m— 

An approach based on logarithmic norms: (compare Lemma fTTT) We fix 
any norm || • ||, for example the L°°-norm: |jx||oo = max, \xi\ (one should chose 
the norm which gives the smallest I ) 

1. [5] = [{f(x,y c ) - f(x,y) \ x € \Wj],y € [W y }}]j. 

2. C=\\[6]\\ 

3. l = nghtU(^([W 2 ],y c ))] 



_ C(e_[ 



-1) 



4. If I ^ 0, then D = 
If / = 0, then D = Ch~ 

5. [A] = [-D,D] n 

Remark. In both cases we compute 

[5] = [{f(x, y c ) - f(x, y) \ x G [W x ], y e 



(59) 



One need to be very careful in the computation of [S] using (|59p , because direct 
interval evaluation of [{/(x, y c ) — f(x,y) | x £ [Wi],y € [Wy]}]/ yields big 
overestimation. Namely, when there is no perturbations at all, i.e. [W y ] = {y c }, 
then [5] = 0. On the other hand if /([W x ]) = [{f(x,y c ) \ x € [Wi]}]/ = [a~,a+} 
then the naive interval computation give [5] = [a - — a + , a + — a - ], so diam [8] — 
2diam/([Wi]) and this can be big because [Wi] is an enclosure of a solution 
during the whole time step. 
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6.4 Rearrangement 

The rearrangement is an essential ingredient in the Lohner algorithm, designed 
to reduce the wrapping effect [Loi ILoll IMoj . We will not discuss this issue here, 
but we will only include necessary formulas (see |Z1] for more details and the 
motivation) . 

Evaluations 2 and 3. In this representation 

[x k ]=x k + [B k ][f k ]. (60) 

In the context of our algorithm in part 3 we obtain 

[x k +i] = x k +i + [B k+ i][r k+1 ]. (61) 

Now we have to take into account equation ([56^1 . We set 

x k+ i = m{x k+1 + [A}) (62) 
[f fe +i] = [f k+1 ] + [B^ 1 ](x k+1 + [A}-x k+1 ). (63) 

Evaluation 4. In this representation 

[x k ]=x k + C k [ro] + [B k ]{f k }. (64) 

In the context of our algorithm in part 3 we obtain 

[x k+ i] = x k +i + C k+ i [r ] + [-Bfc+i] [rk+i] ■ (65) 

Equation (|56J) is taken into account exactly in the same way as in previous 
evaluations, i.e., we use (f6"2"]) and (|6"3"|) . 

7 Rigorous estimates between time steps 

In order to compute the Poincare map for differential inclusion we also need an 
estimate for time t £ [t k ,t k + h k ]. 

Input parameters: 

• h k is a time step 

• [xk] C M™, such that ip(t k , [x ], [yo]) C [x k ] 

• [x k+ x] c R n , such that (p(t k + h k , [x ], [yo]) C [x k+1 ] 

• convex and compact set [Wb] C K n and convex set [W y ] C M m , such that 

<p(.[tk,tk + h k ],[x ],\3fo]) c [W 2 ] (66) 
[y }([t k ,t k+1 ]) c [W v ]. (67) 
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• Vc € [W„] 

• [x k+ i] C R™, such that Tp(h k , [x k ],y c ) C [Sfe+i] 

• [Wi] C R™ compact and convex, such that ^([0, /ifc], [xk],y c ) C [Wi] 
Output: 

We compute C R" such that 

ip(t k + [0,h k ],[x ],[y ]) c 

Algorithm: 

• We compute [-Efc] C R", such that 

^([0,fcfc] ) [a fc ],tf c )c[S fc ] (68) 

using a procedure for an ODE described in |Z1| . This procedure requires 
as input data: h k , [x k ], pSfc+i] an d [Wi]. 

• we compute a set [A] C R™, such that 

<p(tk + h, [xo], [yo]) C Tp(h, [xk],y c ) + [A], for < /i < /i fe . (69) 
Observe that this requires y c , [Wi], [W2] and \W V \. 

• finally we obtain 

v (t k + [0, h k ], [x Q ], [y ])i C [E k ]i = [E k ]i + [A],. (70) 

Slightly better algorithm: 

• if ^ /i([W2], [Wj,])^ then the i-th coordinate is strictly monotone on 
[W2] x [W v ], hence we set 

[E k ]i = hull([x fc ]j, [x k+ i]i) 

• if € fi([W 2 ], [W v ]), then we compute \E k ] C R", such that 

<p([0,h k ],[x k ],y c )c[E k ] (71) 

using a procedure for an ODE described in [Zlj . This procedure requires 
as input data: h k , [x k ], [aJfc+iJand [Wi\. 

We have 

<p(t k + [0, h k ], [x ], [y ])i C [E k ]t = \E k ]i + [A];. (72) 

A drawback of this approach: 

if we have to perform several time steps during which the computed enclosure 
for the trajectory has a nonempty intersection with the section, then A is added 
twice. 
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7.1 Computation of the Poincare map 

If as in |Zlj we assume that the section is given by a(x) =0 then an algorithm 
discussed in Section 5 in |Z1| also applies in the present context. 

8 Some tests, discussion 

8.1 Perturbed harmonic oscillator 

We use the harmonic oscillator to compare two methods: first based on the 
logarithmic norms and the second one that uses component-wise estimates. To 
shorten the notation in this section we call them LN method and CW method 
correspondingly. 

The equations of the perturbed harmonic oscillator are given by 

x' = y + e 1 (73) 

V = -x + e 2 

and we will always use the initial condition given by (1, 0) + [—6, S] 2 . 

In both methods we first find the solution of the unperturbed system and 
then we add the influence of perturbation denoted (following section \§§ by A. 
For this simple system we are able to compute A for both methods by hand. 
Let h denote time step used. 

For LN method we used the euclidian logarithmic norm fj, e because it is optimal 
for this case. Namely, we have 

i = Me(^([W 2 ],»c)=Mef 

Therefore, we obtain A = [-D, D] 2 where 

D = hy/el + e 2 2 . 
For CW method we obtain A = ([— £>i, Z?i], [—D 2: D 2 }), where 

Di = t\ sinh h + £2(cosh h — 1), 
D2 = ei (cosh h — 1) + €2 sinh h. 

Suppose that t\ = €2 '•= e, then LN method is better than CW method if 

V2he < e(sinh h + cosh h - 1) = e(exp(/i) - 1) (75) 

Inequality fTS]) holds for h > 0.657275. As it can be seen in Table Q] results of 
computations agree with this theoretical estimate and the LN method is better 
for h > 0.657275. We were not able to use time steps h > 0.8 because for such 
a big time steps our rough enclosure procedure (the first part of the algorithm) 
fails. 
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timp sfpn 


LN method 


CW method 


h 




^ 1 •> ^ A 


0.799 


0.112996 


0.122332 


0.7 


0.0989949 


0.101375 


0.66 


0.0933381 


0.0934792 


0.658 


0.0930553 


0.0930927 


0.657 


0.0929138 


0.0928997 


0.65 


0.0919239 


0.0915541 


0.5 


0.0707107 


0.0648721 


0.25 


0.0353553 


0.0284025 


0.1 


0.0141421 


0.0105171 


0.01 


0.00141421 


0.00100502 


0.001 


0.000141421 


0.00010005 



Table 1: Perturbed harmonic oscillator ex = e 2 = 0.1: Estimates of perturba- 
tions for various time steps - comparison between LN and CW method 



The situation is quite different, when we perturb only one coordinate. Sup- 
pose that ei = and e 2 = £• Now, for LN method we have 



D = he, 



and for CW method 



h 2 h 4 

Dx = e(cosh/i-l) = e(— + — + ...), 

h? h 5 

D 2 = esmhh = e{h + — + — + ...). 

3! 5! 

From the above formulas it follows that for time steps up to 1.616137 value of 
Dx is smaller than D, but D 2 is always bigger than D. In Table [5] we list values 
of perturbations for LN an CW method for various time steps. Again for time 
steps bigger than 0.8 our implementation could not find rough enclosure. For 
small time steps the ratio jy- is quite big, when the ratio is slightly less than 
one. So overall it is better to use CW method. 

In Table[3]we compare diameters of computed rigorous estimates of solutions 
of (f?3"|) after time T = 2tt for these two methods using various values of h, e and 
S. Again we perturb only second coordinate i.e. ex = 0, e 2 — e. As expected, 
we see that decreasing time steps results in the increase of the accuracy of the 
estimates, but it also increases computational cost. In the second part of the 
table we were changing set sizes and in the third one we were changing the size 
of the perturbation. It can be seen that our algorithm is capable to provide 
estimates even for perturbations much bigger than values of the vector field. 
Observe that with the time steps used in these experiments the CW method is 
better than LN method. The biggest time step h used was approximately equal 
to 0.785. 
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time step 


LN method 


CW method 


h 


D 


D 1 


D 2 


0.8 


0.08 


0.0337435 


0.0888106 


0.5 


0.05 


0.0127626 


0.0521095 


0.25 


0.025 


0.0031413 


0.0252612 


0.1 


0.01 


0.0005004 


0.0100167 


0.01 


0.001 


5.0e-06 


0.0010001 


0.001 


0.0001 


5.002e-08 


0.0001 



Table 2: Perturbed harmonic oscillator e\ = 0,62 = 0.1: Estimates of per- 
turbations for various time steps - comparison between LN and CW method 







number 


size of the set after time T = 2ir 


e 


5 


of steps 


LN method 


CW method 


0.1 


0.01 


8 


1.5789308 


1.2143687 


0.1 


0.01 


100 


1.6220657 


0.8479880 


0.1 


0.01 


1000 


1.6202468 


0.8227680 


0.1 


0.01 


10000 


1.6200250 


0.8202765 


0.1 


0.01 


100000 


1.6200025 


0.8200276 


0.1 





100 


1.5994735 


0.8253958 


0.1 


0.01 


100 


1.6220657 


0.8479880 


0.1 


0.1 


100 


1.8253953 


1.0513176 


0.01 


0.01 


100 


0.1825395 


0.1051317 


0.1 


0.01 


100 


1.6220657 


0.8479880 


1 


0.01 


100 


16.017328 


8.2765505 


10 


0.01 


100 


159.96995 


82.562176 



Table 3: Perturbed harmonic oscillator e± = 0, e 2 = e: Estimates of perturba- 
tions for various values of the parameters - comparison between LN and CW 
method 
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initial set [X] 
perturbations [e] 


(0.0,-10.3, 0.03) + {0} x [-10" 4 , lO" 4 ] 2 

[-io- 4 ,io- 4 ] 3 


P([X}) 
diam P([X]) 




( [-0.3136278,0.3049910] \ T 

[-3.7425421,-3.4205722] 
^ [0.0306989,0.0337781] / 
0.6186189, 0.3219698,0.0030791) 



Table 4: Perturbed Rossler equation: Value of a Poincare map on section = 
{x = 0,x' > 0} 



8.2 Rossler equations 

Rossler equations [RJ are given by 

x' = -{y + z) 

y' = x + 0.2y (76) 
z' = 0.2 + z(x-a), 

where a is a real parameter. In our tests we set a = 5.7 - the 'classical' parameter 
value for which numerical simulation display a strange attractor [R] , 

In our test we focus on computation of a Poincare map, P, on section = 
{x = 0, x' > 0} around a point x = (0.0,-10.3, 0.03). This is a point from the 
attractor (or close to the attractor, which we have found numerically difficult 
in [Z3]). 

In Table [4] we list the results of a computation of Poincare map on section 
for a differential inclusion x' <G f(x) + [e], where f(x) is the vector field 
in Rossler equations (|76[) and [e] = [— 10~ 4 , 10 -4 ] 3 . The initial condition was 
xo + {0} x [— 10~ 4 , 10~ 4 ] 2 . In computations the method based on the component- 
wise estimates and the Lohner algorithm - 4th evaluation was used. 

We see that our algorithm can provide good estimetes even for perturbed 
system and for set of initial data containing numerically difficult points from 
attractor. 

8.3 Kuramoto-Sivashinsky PDE's 

Assuming odd and periodic boundary conditions the Kuramoto-Sivashinsky 
equations can be reduced [ZMj to the following infinite system of ordinary dif- 
ferential equations 

k— 1 oo 

dfc = fc 2 (l - vk 2 )a k - fc^ a n a k - n + 2fcJ^ a n a n+k k= 1,2,3,... (77) 

n— 1 n— 1 

where v > 0. In [Z21IZ4] using the algorithm based on component-wise estimates 
described in this paper to handle the dominant modes and the method of self- 
consistent bounds developed in [ZMj to deal with the tail (the remaining modes) 
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the existence of multiple periodic orbits has been proved for a range for v G 
[0.032,0.127]. Some of these orbits were attracting, while others were unstable 
with one unstable direction. 
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